﻿///*
// ядро NS
//*/
//#include "math.h"
//#include "stdio.h"
//#include "stdlib.h"
//#include "string.h"
//#include "alloc.h"
//#include "conio.h"
//#include "time.h"
//
//#include "nsglobal.h"
//
//ldouble  a[MAXR+1]; /* коэффициенты полиномиального множителя ВЯ */
//
//hpdouble grm;     /* матрица Грама */
//hpdouble p[MAXM]; /* узлы хаотической сетки */
//hpdouble u;       /* значения в этих узлах */
///**/
///*  nsinem  */
//hpdouble mgrk;     /* матрица Грама активного набора */
//hpdouble mv, mu, mc;
//hpint mind;
///**/
//hpdouble xp;
///**/
//hpint vaiw;
///**/
//
//void main()
//{
// int m, n, nfun, r, ti, lit;
// int i, j, nst, ier;
// unsigned int lw;
// double x0, y0, x1, y1, xk, yk, v, vf;
// double eps;
// double ERR;
// int NSINEM_OUT;
// int ij;
// double errp, errs, errm, errx;
// const char *nfin="nsparam.c", *nfout="nsreport.c";
// const char *nfpnt="pnt.dat", *nfmod="mod.dat", *nfres="res.dat",
//            *nfdiff="diff.dat";
// char *ws="                                                                     ";
// FILE *fin, *fout;
// FILE *fpnt, *fmod, *fres, *fdiff;
// clock_t t1, t2;
// double dt;
// double cond, wd;
//
// ier = 0;
// fout = fopen(nfout, "w");
// if( fout == NULL ) {
//   printf("ошибка при открытии файла результатов");
//   getch();
//   exit(0);
// }
//
// fin = fopen(nfin, "r");
// if( fin == NULL ) {
//   fprintf(fout,"ошибка при открытии файла параметров");
//   exit(0);
// }
//
// fpnt = fopen(nfpnt, "w");
// if( fpnt == NULL ) {
//   fprintf(fout,"ошибка при открытии файла узлов");
//   exit(0);
// }
//
// fmod = fopen(nfmod, "w");
// if( fmod == NULL ) {
//   fprintf(fout,"ошибка при открытии файла модели");
//   exit(0);
// }
//
// fres = fopen(nfres, "w");
// if( fres == NULL ) {
//   fprintf(fout,"ошибка при открытии файла сплайна");
//   exit(0);
// }
//
// fdiff = fopen(nfdiff, "w");
// if( fdiff == NULL ) {
//   fprintf(fout,"ошибка при открытии файла отклонений");
//   exit(0);
// }
//
///* r n nfun m eps */
// fgets(ws,strlen(ws)-1,fin);
// fscanf(fin,"%d %d %d %d %lf %d %d %lf %d %d",
//             &r,&n,&nfun,&m,&eps,&nst,&ti,&ERR,&lit,&NSINEM_OUT);
// fprintf(fout,
// "r=%d  n=%d  nfun=%d  m=%d  eps=%le  nst=%d  ti=%d  ERR=%lf  lit=%d  NSINEM_OUT=%d\n",
//  r,n,nfun,m,eps,nst,ti,ERR,lit,NSINEM_OUT);
// switch (ti) {
// case 1  :
//          fprintf(fout,"\nнормальные сплайны : интерполяция\n");
//          break;
// case 2  :
//          fprintf(fout,"\nнормальные сплайны : равн.сглаживание (NSINEM) \n");
//          break;
// case 10 :
//          fprintf(fout,"\nнормальные сплайны : оценка обусловленности матрицы Грама\n");
//          break;
// case 11 :
//          fprintf(fout,"\nобобщенный сплайн (Василенко)\n");
//          break;
// case 12 :
//          fprintf(fout,"\n'иррациональная интерполяция'\n");
//          break;
// default :
//          fprintf(fout,"\nтип интерполяции не определен\n");
// };
//
// lw = m;
// lw = (lw+1)*lw/2;
// if( ti == 11 ) {
//    lw = m + 3;
//    vaiw = (hpint)farmalloc(sizeof(int)*lw);
//    if( vaiw == NULL ) {
//       fprintf(fout,"ошибка при выделении памяти N 1v");
//       exit(0);
//     };
//     lw = m + 3;
//     lw *= lw;
// };
// if( ti == 10 ) {
//   lw = m;
//   lw = lw*lw;
// };
// grm = (hpdouble)farmalloc(sizeof(double)*lw);
// if( grm == NULL ) {
//     fprintf(fout,"ошибка при выделении памяти N 1");
//     exit(0);
// };
// if( ti == 2 ) {
///* NSINEM */
//   lw = m;
//   lw = (lw+1)*lw/2;
//   mgrk = (hpdouble)farmalloc(sizeof(double)*lw);
//   if( mgrk == NULL ) {
//     fprintf(fout,"ошибка при выделении памяти N 1(NSINEM)");
//     exit(0);
//   };
//
//   lw = m;
//   lw *= 4;
//   mu = (hpdouble)farmalloc(sizeof(double)*lw);
//   if( mu == NULL ) {
//     fprintf(fout,"ошибка при выделении памяти N 2(NSINEM)");
//     exit(0);
//   };
//
//   lw = m;
//   lw *= 2;
//   mv = (hpdouble)farmalloc(sizeof(double)*lw);
//   if( mv == NULL ) {
//     fprintf(fout,"ошибка при выделении памяти N 3(NSINEM)");
//     exit(0);
//   };
//
//   lw = m;
//   lw *= 2;
//   mc = (hpdouble)farmalloc(sizeof(double)*lw);
//   if( mc == NULL ) {
//     fprintf(fout,"ошибка при выделении памяти N 4(NSINEM)");
//     exit(0);
//   };
//
//   lw = m;
//   lw *= 3;
//   mind = (hpint)farmalloc(sizeof(double)*lw);
//   if( mind == NULL ) {
//     fprintf(fout,"ошибка при выделении памяти N 5(NSINEM)");
//     exit(0);
//   };
// };
//
// lw = n;
// for(i=0; i<m; i++) {
//   p[i] = (hpdouble)farmalloc(sizeof(double)*lw);
//   if( p[i] == NULL ) {
//     fprintf(fout,"ошибка при выделении памяти N 2");
//     exit(0);
//   }
// }
//
// lw = m;
// if( ti == 11 ) lw = m + 3;
// if( ti == 10 ) lw = 2*m;
// u = (hpdouble)farmalloc(sizeof(double)*lw);
// if( u == NULL ) {
//   fprintf(fout,"ошибка при выделении памяти N 3");
//   exit(0);
// }
//
// lw = m + n;
// xp = (hpdouble)farmalloc(sizeof(double)*lw);
// if( xp == NULL ) {
//   fprintf(fout,"ошибка при выделении памяти N 4");
//   exit(0);
// }
//
// nsfa(r, &(a[0]));
// fprintf(fout,"\n Коэффициенты полиномиального множителя ВЯ \n");
//#ifdef LD
// for(i=0; i<=r; i++) fprintf(fout,"a[%d] = %Le ",i,a[i]);
//#else
// for(i=0; i<=r; i++) fprintf(fout,"a[%d] = %le ",i,a[i]);
//#endif
// fprintf(fout,"\n");
///* расчет оценки наименьшего допустимого значения eps */
// if ( r == 0 ) wd = MEPS;
// else {
//   wd = MEPS*fac(r+1,2*r)/pow(2,r);
//   wd = pow(wd,1./(double)r);
// };
// fprintf(fout,"Оценка наименьшего допустимого значения eps = %le\n",wd);
///**/
///* формирование узлов */
// m = fp(n,m);
///* формирование значений функции в узлах */
// if(ti != 11) fu(m, nfun, u, ERR);
// else vafu(m, nfun, u, ERR);
///**/
// fprintf(fout,"Узлы хаотической сетки \n");
// for(i=0; i<m; i++) {
//   for(j=0; j<n; j++) fprintf(fout,"p[%d][%d]=%le ",i,j,p[i][j]);
//   fprintf(fout," u=%+le \n",u[i]);
// }
// t1 = clock();
// if( ti <= 2 ) nsfgrm(r, n, m, eps, grm, &(a[0]));
// if( ti == 10 ) nsfgrm2(r, n, m, eps, grm, &(a[0]));
// if( ti == 11 ) vafgrm(n, m, grm);
// t2 = clock();
// dt = (t2-t1)/CLK_TCK;
// fprintf(fout,"\nвремя формирования матрицы Грама : %le (с)  %ld %ld\n",
//              dt,t1,t2);
///**/
// t1 = clock();
// if (ti == 1 ) ier = holds(m, grm);
// if( ier != 0 ) {
//   fprintf(fout,"матрица Грама вырождена\n");
//   exit(0);
// }
// if(ti == 1) hols(m, grm, u);
// if (ti == 2) {
//    int li;
//    int mkp, m1, m2;
//    FILE *mfl;
//    const char *nmfl="nsinem.c00";
//
//    mfl = fopen(nmfl, "w");
//    if( mfl == NULL ) {
//      printf("ошибка при открытии файла результатов NSINEM");
//      getch();
//      exit(0);
//    };
//
//    m1 = 0;
//    m2 = m;
//
//    for(i=1;i<=m1;i++) {
//     	mind[i-1] = i;
//	    mc[i-1] = u[i-1];
//    };
//    if(ERR <= 0.) {
//      ERR = 0.1;
//      fprintf(fout,"NSINEM : ERR = %lf\n",ERR);
//    };
//    for(i=1;i<=m2;i++) {
//	    mc[m1+i-1] =  u[m1+i-1] + ERR;
//	    mc[m+i-1]  = -u[m1+i-1] + ERR;
//  	};
//
//    lw = (int)(m+1)*m/2;
//    for(li=0; li<lw; li++) mgrk[li] = grm[li];
//    ier = holds(m, mgrk);
//    if( ier != 0 ) {
//      fprintf(fout,"матрица Грама вырождена\n");
//      exit(0);
//    };
//    for(i=0; i<m; i++) mv[i] = u[i];
//    hols(m, mgrk, mv);
//
//    if(!NSINEM_OUT) mfl = NULL;
//    ier = nsinem(m1, m2, &mkp, mind, grm, mgrk, mv, mu, mc, &lit, mfl);
//
//    fclose(mfl);
//    fprintf(fout,"NSINEM : mkp = %d\n",mkp);
//    fprintf(fout,"NSINEM : lit = %d\n",lit);
//    fprintf(fout,"NSINEM : ier = %d\n",ier);
//    if( ier == 0 ) fprintf(fout,"NSINEM : O'K\n");
//    if( ier == 2 ) fprintf(fout,"NSINEM : Лимит итераций\n");
//    if( ier == 3 ) {
//      fprintf(fout,"NSINEM : Ошибка в данных\n");
//      exit(0);
//    };
//    if( ier == 1 ) {
//      fprintf(fout,"NSINEM : матрица Грама вырождена\n");
//      exit(0);
//    };
//    for(i=0; i<m; i++) u[i] = mv[i];
// };
// if (ti == 11) {
//   cond = decomp(m+3, grm, vaiw, xp);
//   fprintf(fout,"\nDECOMP : cond = %le\n",cond);
//   solve(m+3, grm, u, vaiw);
// };
// if (ti == 10) {
//    int irang;
//    double tola, tolb;
//    tola = tolb = 0.;
//    irang = 0;
//
//    ier = evdd(m, grm, u);
//
//    if( ier == 0 ) ier = evdp(0, m, u, tola, tolb, cond, irang);
//    if( ier == 0 ) {
///* получим (точную) оценку обусловленности для grm */
//       cond = 0.;
//       j = m-1;
//       while( j >= 0 && u[j] <= 0. ) j--;
//       if( j > 0 && u[j] > 0. ) cond = u[0]/u[j];
//       ier = evdp(1, m, u, tola, tolb, cond, irang);
//    };
//    fprintf(fout,"\nEVDD : ier = %d   cond = %le\n",ier,cond);
// };
// t2 = clock();
// dt = (t2-t1)/CLK_TCK;
// fprintf(fout,"\nвремя определения коэф.сплайна : %le (с)  %ld %ld\n",
//               dt,t1,t2);
///**/
// if(ti == 10) goto lkex;
//
// fprintf(fout,"\nЗначения в узлах хаотической сетки \n");
// if(ti != 11) fu(m, nfun, xp, ERR);
// else vafu(m, nfun, xp, ERR);
// errp = 0.;
// fprintf(fpnt,"%d %d\n",m,n);
// for(i=0; i<m; i++) {
//   for(j=0; j<n; j++) {
//    fprintf(fout,"p[%d][%d]=%le ",i,j,p[i][j]);
//    fprintf(fpnt," %le ",p[i][j]);
//   };
//   switch (ti) {
//   case 1  :
//   case 2  : v = nsv(r, m, n, eps, p[i], u, &(a[0])); break;
//   case 11 : v = vav(m, n, p[i], u); break;
//   case 12 : v = irv(m, n, p[i], u); break;
//   default : v = 0.;
//   };
//   errp += fabs(xp[i] - v);
//   fprintf(fout,"функ.=%+le  спл.=%+le  функ.-спл.=%+le \n",xp[i],v,xp[i]-v);
//   fprintf(fpnt," %+le \n",xp[i]);
// }
// fprintf(fout,"\nсуммарное нарушение в узлах хаотической сетки : %le\n",errp);
//
// fprintf(fout,"\nЗначения функции\n");
//
// fprintf(fmod," %d %d\n",nst+1,nst+1);
// fprintf(fres," %d %d\n",nst+1,nst+1);
// fprintf(fdiff," %d %d\n",nst+1,nst+1);
//
// x0 = y0 = 0.;
// x1 = y1 = 1.;
// fprintf(fout,"\n%9s","         ");
// for(j=0; j<nst+1; j++){
//    xk = x0 + j*(x1 - x0)/nst;
//    fprintf(fout," %lf ",xk);
//    fprintf(fmod," %lf ",xk);
// };
// fprintf(fout,"\n");
// fprintf(fmod,"\n");
// for(i=0; i<nst+1; i++) {
//   yk = y1 - i*(y1 - y0)/nst;
//   fprintf(fout,"%lf ",yk);
//   fprintf(fmod,"%lf ",yk);
//   for(j=0; j<nst+1; j++){
//      xk = x0 + j*(x1 - x0)/nst;
//      switch (nfun) {
//      case 1 : vf = tf1(xk, yk); break;
//      case 2 : vf = tf2(xk, yk); break;
//      case 3 : vf = tf3(xk, yk); break;
//      default : vf = 0.;
//      };
//      fprintf(fout,"%+lf ",vf);
//      fprintf(fmod,"%+lf ",vf);
//  };
//  fprintf(fout,"\n");
//  fprintf(fmod,"\n");
// };
//
// errs = 0.;
// errm = 0.;
// errx = 0.;
// ij = 0;
// fprintf(fout,"\nЗначения сплайна\n");
// x0 = y0 = 0.;
// x1 = y1 = 1.;
// fprintf(fout,"\n%9s","         ");
// for(j=0; j<nst+1; j++){
//    xk = x0 + j*(x1 - x0)/nst;
//    fprintf(fout," %lf ",xk);
//    fprintf(fres," %lf ",xk);
//    fprintf(fdiff," %lf ",xk);
// };
// fprintf(fout,"\n");
// fprintf(fres,"\n");
// fprintf(fdiff,"\n");
//
// t1 = clock();
//
// for(i=0; i<nst+1; i++) {
//   yk = y1 - i*(y1 - y0)/nst;
//   fprintf(fout,"%lf ",yk);
//   fprintf(fres,"%lf ",yk);
//   fprintf(fdiff,"%lf ",yk);
//   for(j=0; j<nst+1; j++){
//      xk = x0 + j*(x1 - x0)/nst;
//      xp[0] = xk;
//      xp[1] = yk;
//      switch (ti) {
//      case 1  :
//      case 2  : v = nsv(r, m, n, eps, xp, u, &(a[0])); break;
//      case 11 : v = vav(m, n, xp, u); break;
//      case 12 : v = irv(m, n, xp, u); break;
//      default : v = 0.;
//      };
//      switch (nfun) {
//      case 1 : vf = tf1(xk, yk); break;
//      case 2 : vf = tf2(xk, yk); break;
//      case 3 : vf = tf3(xk, yk); break;
//      default : vf = 0.;
//      };
//      errs += (v - vf)*(v -vf);
//      errm += fabs(v - vf);
//      if ( fabs(v - vf) > errx ) errx = fabs(v - vf);
//      ij++;
//      fprintf(fout,"%+lf ",v);
//      fprintf(fres,"%+lf ",v);
//      fprintf(fdiff,"%+lf ",v-vf);
//  };
//  fprintf(fout,"\n");
//  fprintf(fres,"\n");
//  fprintf(fdiff,"\n");
// };
//
// t2 = clock();
// dt = (t2-t1)/CLK_TCK;
// fprintf(fout,"\nвремя расчета зн.сплайна в узлах рег.сетки : %le (с)  %ld %ld\n",
//              dt,t1,t2);
//
// errs = sqrt(errs/ij);
// errm = errm/ij;
// fprintf(fout,"\nср.кв. отклонение в узлах регулярной сетки : %le\n",errs);
// fprintf(fout,"\nср.аб. отклонение в узлах регулярной сетки : %le\n",errm);
// fprintf(fout,"\nма.аб. отклонение в узлах регулярной сетки : %le\n",errx);
//
//lkex :
//
// farfree((void far*)grm);
// for(i=0; i<m; i++) farfree((void far*)p[i]);
// farfree((void far*)u);
// farfree((void far*)xp);
// if( ti == 11 ) farfree((void far*)vaiw);
// if( ti == 2 ) {
//  farfree((void far*)mgrk);
//  farfree((void far*)mu);
//  farfree((void far*)mv);
//  farfree((void far*)mc);
//  farfree((void far*)mind);
// };
// fprintf(fout,"\n ier = %d\n",ier);
//
// fclose(fin);
// fclose(fpnt);
// fclose(fmod);
// fclose(fres);
// fclose(fdiff);
// fclose(fout);
//
// return;
//}
//
//double tf1(double x, double y)
///* тестовая функция N1 (Василенко) */
//{
// double f;
//
// f = 3*x*x*sin(5*y+x) + 4*(y-x)*cos(3*x+2*y);
// return f;
//}
//
//double tf2(double x, double y)
///* тестовая функция N2 (Пикалов) */
//{
// double f;
// double s, SMIN;
// double a[3], b[3], c[3], x0[3], y0[3];
// int i;
//
// SMIN = -100.;
//
// a[0] = a[2] = 9.;
// a[1] = 4.5;
// b[0] = b[1] = b[2] = 9.;
// c[0] = 5.;
// c[1] = 0.;
// c[2] = 2.5;
// x0[0] = x0[1] = x0[2] = 0.5;
// y0[0] = 0.25;
// y0[1] = 0.;
// y0[2] = 0.75;
//
// f = 0.;
// for(i=0; i<3; i++) {
//   s = -a[i]*a[i]*(x-x0[i])*(x-x0[i]) - b[i]*b[i]*(y-y0[i])*(y-y0[i]);
//   if ( s < SMIN ) s = SMIN;
//   f += c[i]*exp(s);
// }
// return f;
//}
//
//double tf3(double x, double y)
///* тестовая функция N3 */
//{
// double f, r, v, x0, y0;
//
// x0 = y0 = 0.5;
// r = 0.25 + 10.*MEPS;
// v = 1.0 + 10.*MEPS;
// f = 0.;
// if ( (x-x0)*(x-x0)+(y-y0)*(y-y0) < r*r ) f = v;
// return f;
//}
//
//void fu(int m, int nf, hpdouble u, double ERR)
///* формирование значений тестовой функции в узлах хаотической сетки */
//{
// int i, ini;
// double w;
//
// ini=0;
// srand(ini);
// for(i=0; i<m; i++) {
//      switch (nf) {
//      case 1 : u[i] = tf1( p[i][0], p[i][1]); break;
//      case 2 : u[i] = tf2( p[i][0], p[i][1]); break;
//      case 3 : u[i] = tf3( p[i][0], p[i][1]); break;
//      default : u[i] = 0.;
//      };
//      w = (double)rand()/(double)RAND_MAX;
//      u[i] += ( 2.*ERR*(w - 0.5) );
// };
// return;
//}
//
//int fp(int n, int m)
///* формирование узлов хаотической сетки */
//{
// int i,j, ini;
//
// if (m>MAXM-1) m=MAXM-1;
//
// ini=666;
// srand(ini);
// for(i=0; i<m; i++)
//   for(j=0; j<n; j++) p[i][j] = (double)rand()/(double)RAND_MAX;
// return m;
//}
//
//double nsggrm(int r, int n, double eps, int i, int j, ldouble *a)
///* формирование элемента матрицы Грама, функционалы - значения функции
//   в узлах */
//{
// double t, v;
//
// t = nsft(n, eps, p[i], p[j]);
// v = nsrk(r, t, a);
// return v;
//}
//
//void nsfgrm(int r, int n, int m, double eps, hpdouble grm, ldouble *a)
///* формирование матрицы Грама, функционалы - значения функции в узлах */
//{
// int i, j, ij;
//
// ij = 0;
// for(i=0; i<m; i++) {
//   for(j=i; j<m; j++) {
//    grm[ij] = nsggrm(r, n, eps, i, j, a);
///*
//    t = nsft(n, eps, p[i], p[j]);
//         fprintf(fout,"i =%d  j =%d  ij =%d  t =%le \n",i,j,ij,t);
//    grm[ij] = nsrk(r, t, a);
//*/
//    ij++;
//   }
// }
// return;
//}
//
//void nsfgrm2(int r, int n, int m, double eps, hpdouble grm, ldouble *a)
///* формирование полной (квадратной) матрицы Грама,
//   функционалы - значения функции в узлах */
//{
//#define A(I_,J_) (grm+(int)(I_)*(m)+(J_))
// int i, j;
//
// for(i=0; i<m; i++)
//   for(j=i; j<m; j++) {
//    *A(i,j) = nsggrm(r, n, eps, i, j, a);
//    *A(j,i) = *A(i,j);
//   }
//
// return;
//#undef A
//}
//
//
//double phi0(int n, hpdouble x)
///* формирование значения пробной функции в x */
//{
// double f;
//
// f = 0.;
// return f;
//}
//
//double nsv(int r, int m, int n, double eps, hpdouble x, hpdouble u, ldouble* a)
///* вычисление значения сплайна в x */
//{
// double t, v;
// int i;
//
// v = phi0(n,x);
// for(i=0; i<m; i++) {
//    t = nsft(n, eps, p[i], x);
//    v += u[i]*nsrk(r, t, a);
// };
// return v;
//}
//
///**/
///**/
//
//double fac(int k1, int k2)
///* вычисление произведение k1*(k1+1)...(k2-1)k2 */
//{
// int i;
// double f;
//
// f = 1.;
// if( k2 <= 1 ) return f;
// for(i=k1; i<=k2; i++) f *= i;
// return f;
//}
//
//void nsfa(int r, ldouble* a)
///* формирование коэффициентов полиномиального множителя ВЯ;
//   минимизация вычислительных погрешностей */
//{
// int k, i;
// double s1, s2;
//
// a[r] = 1.;
// if (r == 0) return;
// a[r-1] = 1.;
// if (r == 1) return;
//
// for(k=0; k<=r-2; k++) {
//   s1 = 1.;
//   for(i=1; i<=r-k; i++) s1 *= ((ldouble)2./i);
//   s2 = 1.;
//   for(i=1; i<=r; i++) s2 *= ((ldouble)(k+i)/(r+i));
//   a[k] = s1 * s2;
// }
//
// return;
//}
//
//void nsfa_old(int r, ldouble* a)
///* формирование коэффициентов полиномиального множителя ВЯ */
//{
// int k;
// double s1;
//
// a[0] = 1.;
// if (r == 0) return;
// s1 = 1.;
// for(k=1; k<=r; k++) {
//   s1 *= 2;
//   a[k] = fac(k+1,r+k)/(s1*fac(1,r-k));
// }
// a[0] /= a[r];
// for(k=1; k<=r; k++) a[k] /= a[r];
// return;
//}
//
//double nsp(int r, double t, ldouble* a)
///* вычисление значения полиномиального множителя ВЯ по схеме Горнера */
//{
// int i;
// double s;
//
// s = a[0];
// for(i=1; i<=r; i++) s = s*t + a[i];
// return s;
//}
//
//double nsrk(int r, double t, ldouble* a)
///* вычисление значения ВЯ в t */
//{
// double s;
//
// s = exp(-t)*nsp(r,t,a);
// return s;
//}
//
//double nsft(int n, double eps, hpdouble eta, hpdouble x)
///* вычисление значения t */
//{
// int i;
// double s,w;
//
// s = 0.;
// for(i=0; i<n; i++) {
//  /*
//  fprintf(fout,"i =%d  eta[i] =%lf  x[i]=%le \n",i, eta[i], x[i]);
//  */
//    w = eta[i] - x[i];
//    s += (w*w);
// };
// s = eps*sqrt(s);
// return s;
//}
//
///* ir */
//
//double irft(int n, hpdouble eta, hpdouble x)
///* вычисление значения t */
//{
// int i;
// double s, w, t;
//
// s = 0.;
// for(i=0; i<n; i++) {
//    w = eta[i] - x[i];
//    s += (w*w);
// };
// t = 1.e+20;
// if (s != 0.) t = (1./s)*(1./s);
// return t;
//}
//
//double irv(int m, int n, hpdouble x, hpdouble u)
///* вычисление значения "иррационального интерполянта" в x */
//{
// double s1, s2, w, val;
// int i;
//
// s1 = 0.;
// s2 = 0.;
// for(i=0; i<m; i++) {
//    w = irft(n, p[i], x);
//    s1 += (u[i]*w);
//    s2 += w;
// };
// val = s1/s2;
// return val;
//}
//
///* va */
//double vaft(int n, hpdouble eta, hpdouble x)
///* вычисление значения t */
//{
// int i;
// double s, w;
//
// s = 0.;
// for(i=0; i<n; i++) {
//    w = eta[i] - x[i];
//    s += (w*w);
// };
// s = sqrt(s);
// return s;
//}
//
//void vafgrm(int n, int m, hpdouble a)
///* формирование матрицы системы для обобщенного сплайна Василенко */
//{
//#define A(I_,J_) (a+(int)(I_)*((m+3))+(J_))
// int i, j;
// double t, w;
//
// for(i=1; i<=m; i++) {
//   *A(i-1,i-1) = 0.;
//   *A(i-1,m+1-1) = 1.;
//   *A(i-1,m+2-1) = p[i-1][0];
//   *A(i-1,m+3-1) = p[i-1][1];
//   for(j=i+1; j<=m; j++) {
//    t = vaft(n, p[i-1], p[j-1]);
//    if (t<=MEPS) w = 0.;
//    else w = t*t*log(t);
//    *A(i-1,j-1) = w;
//   }
// }
//
// for(i=m+1; i<=m+3; i++)
//   for(j=i; j<=m+3; j++) *A(i-1,j-1) = 0.;
//
// for(i=1; i<=m+3; i++)
//   for(j=1; j<i; j++) *A(i-1,j-1) = *A(j-1,i-1);
//
// return;
//#undef A
//}
//
//void vafu(int m, int nf, hpdouble u, double ERR)
///* форм. значений тест.функции в узлах хаот.сетки сплайн Василенко*/
//{
// int i, ini;
// double w;
//
// ini=0;
// srand(ini);
// for(i=0; i<m; i++) {
//      switch (nf) {
//      case 1 : u[i] = tf1( p[i][0], p[i][1]); break;
//      case 2 : u[i] = tf2( p[i][0], p[i][1]); break;
//      case 3 : u[i] = tf3( p[i][0], p[i][1]); break;
//      default : u[i] = 0.;
//      };
//      w = (double)rand()/(double)RAND_MAX;
//      u[i] += ( 2.*ERR*(w - 0.5) );
// };
// u[m+1-1] =0.;
// u[m+2-1] =0.;
// u[m+3-1] =0.;
// return;
//}
//
//double vav(int m, int n, hpdouble x, hpdouble u)
///* вычисление значения сплайна Василенко в x */
//{
// double t, w, v, st;
// int i;
//
// v = 0.;
// for(i=0; i<m; i++) {
//    t = vaft(n, p[i], x);
//    st = t*t;
//    if (t<=MEPS) w = 0.;
//    else w = 0.5*st*log(st);
//    v += u[i]*w;
// };
// v += u[m+1-1]*1. + u[m+2-1]*x[0] + u[m+3-1]*x[1];
// return v;
//}
//
// 
